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ABSTRACT 

We present an updated mass model for M31 that makes use of a Spitzer 3.6 /im 
image, a mass-to-light ratio gradient based on the galaxy's B ~ R colour profile, and 
observed rotation curve data from a variety of sources. We examine cases where the 
dark matter follows a pure NFW profile and where an initial NFW halo contracts 
adiabatically in response to the formation of the galaxy. We find that both of these 
scenarios can produce a reasonable fit to the observed rotation curve data. However, a 
pure NFW model requires a concentration c v i r = 51 that is well outside the range pre- 
dicted in ACDM cosmology and is therefore disfavoured. An adiabatically contracted 
NFW halo favors an initial concentration c v - 1T = 20 and virial mass 8.2 x 1O 11 M0, and 
this is in line with the cosmological expectations for a galaxy of the size of M31. The 
best-fit mass is consistent with published estimates from Andromeda Stream kinemat- 
ics, satellite galaxy radial velocities, and planetary nebulae studies. Finally, using the 
known linear correlation between rotation curve shear and spiral arm pitch angle, we 
show that the stellar spiral arm pitch angle of M31 (which cannot be deduced from 
imaging data due to the galaxy's inclination) is P = 24?7 ± 4?4. 
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1 INTRODUCTION 

Modelling the mass distribution of M31 is a classical prob- 
lem that has seen many past iterations (e.g., Einasto 1972; 
Kent, Huchra & Stauffer 1989; Klypin, Zhoa & Somerville 
2002; Widrow, Perrett & Suyu 2003; Geehan et al. 2006). 
Recently, a Spitzer 3.6 /im mosaic image of M31 has become 
available (Barmby et al. 2006). This allows for a more ac- 
curate determination of the baryonic mass distribution in 
M31 than has been determined previously. Accurate mass 
modelling of M31 is an important tool for both exploring 
our ideas of cosmology, and for understanding the dynamics 
of newly discovered Local Group dwarf galaxies (e.g., Ma- 
jewski et al. 2007; Chapman et al. 2007) . 

Because M31 can be studied in exquisite detail, it pro- 
vides a crucial testing ground for our ideas of galaxy for- 
mation (Kent 1989; Evans & Wilkinson 2000). A problem 
of particular relevance for Cosmological constant plus Cold 
Dark Matter (ACDM) cosmology is the Tully-Fisher zero- 
point problem, which refers to the fact that standard mod- 
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els cannot reproduce the relation between galaxy luminosity 
and circular velocity (Tully & Fisher 1977) without over- 
producing the number density of galaxies at fixed luminos- 
ity (e.g. Gonzalez et al. 2000; Cole et al. 2000; Benson et 
al. 2003; Yang, Mo & van den Bosch 2003). Another prob- 
lem of related concern is the cusp/concentration problem - 
namely that dark-matter dominated galaxy rotation curves 
seem to rise more slowly than predicted in ACDM (Moore 
1994; Flores & Primack 1994; Moore et al. 1999; van den 
Bosch & Swaters 2001; Blais-Ouellette, Amram & Carignan 
2001; Alam, Bullock & Weinberg 2002; Swaters et al. 2003; 
Simon et al. 2005; Dutton et al. 2005; Kuzio de Naray et al. 
2006). It remains to be seen whether these problems point 
to some unaccounted for process in galaxy formation (e.g. 
Dutton et al. 2007) or to some new physics of cosmological 
relevance (Kaplinghat, Knox & Turner 2000; Zentner & Bul- 
lock 2002; Kaplinghat 2005; Cembranos et al. 2005; Strigari, 
Kaplinghat & Bullock 2007; Gnedin et al. 2007). 

Recently, Dutton et al. (2007) and Gnedin et al. (2007) 
revisited the Tully-Fisher problem using two large well- 
defined samples of disc-dominated (late-type) galaxies. They 
both took as a starting point the "standard" model of disc 
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formation, which assumes that initial Navarro, Frenk & 
White (1997; hereafter NFW) dark haloes respond to disc 
formation via adiabatic contraction (AC) (Blumenthal et al. 
1986; Gnedin et al. 2004; Sellwood & McGaugh 2005; Choi 
et al. 2006). Both groups conclude that the Tully-Fisher zero 
point cannot be explained if initial NFW haloes have concen- 
trations as high as those expected for ACDM with as — 0.9 
(with c ~ 12 for M31-size haloes, e.g. Bullock et al. 2001a, 
Maccio et al. 2006). Both sets of authors agree that the prob- 
lem could be alleviated if AC did not operate (e.g. Somerville 
& Primack 1999) and Dutton et al. (2007) go on to advo- 
cate a model where disc formation induces an expansion in 
the underlying halo density structure. In contrast, Gnedin 
et al. (2007) argue that AC is a fundamental prediction in 
galaxy formation and instead suggest that the cosmology 
be changed to favor lower halo concentrations (e.g. Zentner 
& Bullock 2002) or that the IMF is lighter (with a stellar 
mass-to-light ratio M*/L lower by 0.15 dex) than the stan- 
dard Kroupa (2001) assumption. 

Given the fundamental issues at hand, the question of 
whether AC occurs in nature is of significant interest. The- 
ory certainly favours the idea that haloes contract. Indeed, 
halo contraction must occur when the infall of baryons is 
smooth and adiabatic (e.g. Blumenthal et al. 1986; Ryden 
& Gunn 1987; Sellwood & McGaugh 2005) and simulations 
suggest that dark haloes will contract even when galaxies 
or galaxy clusters form quickly from an irregular collapse 
(Barnes 1987; Flores et al. 1993; Jesseit, Naab & Burkert 
2002; Gnedin et al. 2004; Sellwood & McGaugh 2005; Choi 
et al. 2006; Weinberg et al. 2008). On the other hand, there 
is very little observational evidence that halo contraction 
actually occurs in nature. For example, Zappacosta et al. 
(2006) performed a detailed XMM study of the radio-quiet 
galaxy cluster Abell 2589 and conclude that an NFW halo 
+ AC model cannot explain the data but that a pure NFW 
halo provides a remarkable fit down to ~ 1% of the halo's 
virial radius. Similarly, an investigation of seven elliptical 
galaxies with Chandra by the same group (Humphrey et 
al. 2006) finds that AC degrades the mass profile fits sig- 
nificantly unless strong deviations from a Kroupa IMF are 
allowed. Also, Kassin, de Jong & Pogge (2006a) and Kassin, 
de Jong & Weiner (2006b) found that the rotation curves 
of 32 out of 34 of the bright spiral galaxies they study are 
better fit without adiabatic contraction. 

Klypin et al. (2002) have shown that the observed rota- 
tion curve of M31 is best explained with a halo model that 
includes AC. However, other recent M31 mass modelling pa- 
pers (e.g., Widrow et al. 2003; Geehan et al. 2006) have not 
directly addressed the question of AC. Furthermore, M31 is 
essentially the only well-studied galaxy that has been shown 
to require AC (Klypin et al. 2002). Most other dynamical 
studies of galaxies tend to disfavour AC (e.g., Kassin et al. 
2006a, b). Given the availability of a Spitzer 3.6 /im image, 
it is now appropriate to revisit this question and test the 
conclusion that M31 requires AC. In what follows we show 
that the rotation curve of M31 can only be explained with- 
out AC if an unusually high NFW concentration, c v i r ~ 50, 
is adopted. Such a concentration is well outside the range 
expected in ACDM haloes, and therefore this model is dis- 
favoured. Moreover, the best-fit AC model has a slightly 
high NFW concentration, c v i r = 20, for typical, as = 0.9, 
ACDM haloes (Bullock et al. 2001a, Maccio et al. 2006), 



but it is within 1.5a of the expected value for an M31 size 
halo. We speculate on the implications of these results in 
the conclusion section. 

Our approach is unique compared to other recently pub- 
lished mass profiles derived from modelling the rotation 
curve of M31. While our method bears most similarity to 
that adopted by Klypin et al. (2002), we use an updated 
baryonic mass component, which is derived using a 3.6 /im 
Spitzer image. Geehan et al. (2006) used a pure NFW halo 
profile, and did not consider AC. Tempel et al. (2007) also 
use the 3.6 /im Spitzer image, yet once again, AC is not con- 
sidered. Lastly, Widrow et al. (2003) used a different func- 
tional form for the halo density profile (see Evans 1993). 
The Evans profile is based upon an isothermal distribution 
function. Unlike the NFW profile, the Evans profile was not 
derived by fitting models to dark matter haloes produced in 
ACDM simulations. Widrow et al. (2003) then go on to fit an 
NFW model to their best-fitting Evans (1993) model. How- 
ever, we have chosen to use the cosmologically motivated 
NFW halo density profile model directly. 

Throughout this paper we assume a flat ACDM cos- 
mology with Q m = 0.27 and a Hubble constant of Ho = 75 
km s _1 Mpc~ . We relate virial masses (M v i r ) and radii 
(R v ir) assuming a virial over-density relative to average of 
A vir = 34fl At times we quote baryon fractions relative to 
the universal value with /j, = f2b/f2 m = 0.16. We adopt a 
distance of 784 kpc to M31 (Holland 1998). As a result, an 
angular distance of 1' is equivalent to a distance of 228 pc. 



2 DATA 

We make use of the Spitzer/IRAC 3.6 fj,m image of M31 
(Barmby et al. 2006). The IRAC observations of M31 were 
taken as part of Spitzer General Observer program 3126 in 
2005 August. Fifteen Astronomical Observation Requests 
(AORs) were used to map a region approximately 3? 7 x 1?6 
on the sky. The central 1? 6 x 0? 4 was covered by three AORs, 
each having two 12 s frames per position. The outer regions 
were covered by two AORs, each with two dithered 30 s 
frames per position. The pixel scale of IRAC at 3.6 /im is 
1.221 arcsec, with a point spread function (PSF) that has 
a FWHM ~1.2 arcsec. The 3.6 ^m filter is 1.2 ^m broad, 
spanning the wavelengths 2.965-4.165 fj,m. For a description 
of the IRAC instrument see Fazio et al. (2004). 

The Spitzer data are preferred over 2MASS A" s -band 
data for a few reasons. The regular 2MASS iG-band image is 
simply not deep enough to obtain an accurate bulge/disc de- 
composition. Furthermore, this same 2MASS image suffers 
from over subtraction of the sky by about 0.5 magnitudes 
(T. Jarrett, private communication). A 6x 2MASS AT s -band 
image is also available (where the data were observed with 
an integration time 6 times longer than the regular 2MASS 
survey), but it is only about 1 magnitude deeper and also 
hampered by uncertain sky background. The Spitzer image 
has a larger field of view, better depth, and a more accurate 
determination of the sky level. 

We also use the B — R colour profile of Walterbos & 
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Kennicutt (1987), in order to determine the stellar mass-to- 
light ratio as a function of radius. We have corrected this 
B — R profile has for Galactic foreground extinction, which 
is E(B - V) = 0.062 mags in the direction of M31 (Schlegel, 
Finkbeiner & Davis 1998). This is converted to a correction 
for the B—R colour by using the transformation, E(B—R) = 
IME(B ~ V) from Schlegel et al. (1998), i.e. E(B - R) = 
0.102 for M31. As a comparison the amount of extinction at 
3.6 fim is 1.5 per cent of that at B, and about 34 per cent 
of that in the K„ band. 

We adopt the rotation curve data from several sources. 
In one case (Ml) we adopt the Ha rotation data out to 25 
kpc from Rubin & Ford (19700 

and extend the rotation 
curve to 35 kpc using HI data from Carignan et al. (2006) . 
In another case (M2) we use the CO rotational velocities 
from Loinard, Allen & Lequeux (1995) and the HI from 
Brinks & Burton (1984) to construct an observed rotation 
curve out to a 30 kpc radius. This was the rotation curve 
adopted by Klypin et al. (2002) in their M31 model. We 
adopt Ml as our fiducial case here because Rubin & Ford 
published errors on their observed velocities. However, as we 
show below, both cases yield consistent results. In both cases 
we take into account the presence of a 1O 8 M0 supermassive 
black hole (Bender et al. 2005) and the HI gas distribution 
(Carignan et al. 2006). 



3 MASS MODELLING OF M31 
3.1 The baryonic contribution 

Our goal is to determine the best possible mass model for 
M31. We perform a bulge-disc decomposition in order to 
estimate the baryonic contribution to the rotation curve. 
We then determine several possible mass models. The best 
fit model is determined by minimizing the reduced-^ 2 in a 
fit to the observed rotation curve. 

We first extract the surface brightness profile of M31 
using the Spitzer 3.6 /im image and the IRAF ellipse rou- 
tine, which fits ellipses to an image using an iterative method 
described by Jedrzejewski (1987). In order to mask out fore- 
ground stars and satellite galaxies SExtractor (Bertin & 
Arnouts 1996) was used. An inclination correction was then 
applied to the surface brightness profile (see de Jong 1996; 
Seigar & James 1998a) as follows 



Radius (kpc) 
5 10 15 



2.5Clog(£) 



(1) 



where fii is the surface brightness when viewed at some in- 
clination, i, fi is the corrected surface brightness, a is the 
major axis, b is the minor axis and C is a factor dependent 
on whether the galaxy is optically thick or thin; if C = 1 
then the galaxy is optically thin; if C — then the galaxy is 
optically thick. Since we are using infrared data, the correc- 
tion for inclination is small, and in many cases completely 
ignored (see e.g. Seigar & James 1998a; de Jong 1996). How- 
ever, Graham (2001) showed that C = 0.91 is a good typical 



2 Figure 9 of Rubin & Ford (1970) also includes rotation velocities 
calculated from a narrow [Nil] A6583 emission feature, for the 
region within 3 kpc. However, these velocities are not tabulated 
in their paper, and so we do not include them here. 
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Figure 1. The Spitzer 3.6 /un surface brightness profile of M31 
(without correction for inclination) with decomposition into bulge 
and disc components. The bulge has been fitted with a Sersic 
model (short dashed line) and the disc has been fitted with an 
exponential model (long dashed line). 



value to use for the near-infrared if s -band. If one adopts a 
simple reddening law, where extinction falls as the square of 
wavelength, then the value of C = 0.97 is appropriate at 3.6 
/jm, and we adopt this value here. 

The resulting surface brightness profile (Figure 1) 
reaches a surface brightness of ps.s ~ 21.2 mag arcsec -2 at 
a radius of 1?5. This is very similar to the surface brightness 
found by Barmby et al. (2006). Furthermore, the PNe counts 
of Merrett et al. (2006) predict an R band surface bright- 
ness of ~ 23.8 mag arcsec -2 . At these radii the metallicity of 
M31 is QAZq (Brewer, Richer & Crabtree 1995; Worthey et 
al. 2005). Using the Maraston (2005) population synthesis 
codes, one can see that the expected R— [3.6] color at this 
radius is R— [3.6] ~2. 5, very similar to our measured surface 
brightness. 

From this surface brightness profile, we perform a one- 
dimensional bulge-disc decomposition, which employs the 
Sersic model for the bulge component and an exponential 
law for the disc component (e.g. Andredakis et al. 1995; 
Seigar & James 1998a; Khosroshahi, Wadadekar & Kem- 
bhavi 2000; D'Onofrio 2001; Graham 2001b; Mollenhoff & 
Heidt 2001; see also Graham & Driver 2005 for a review). 
The Sersic (1963, 1968) R}' n model is most commonly ex- 
pressed as a surface brightness profile, such that 



H{R) = /i e exp 



\Re 



1/n 



- 1 



(2) 



where [i e is the surface brightness at the effective radius R e 
that encloses half of the total light from the model (Ciotti 
1991; Caon, Capaccioli & D'Onofrio 1993). The constant b n 
is defined in terms of the parameter n, which describes the 
overall shape of the light profile. When n = 4 the Sersic 
model is equivalent to a de Vaucouleurs (1948, 1959) R 1 ^ 4 
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Figure 2. K— [3.6] colour as a function of B — R colour for Bruzual & Chariot (2003) population synthesis models {Left Panel) and 
Maraston (2005) population synthesis models {Right Panel). In both cases, dashed lines represent a constant age stellar population and 
solid lines represent constant metallicity. The dotted lines represent the typical B — R colour of the inner disc {B — R = 1.79) and the 
outer disc {B - R = 1.40) from Walterbos & Kennicutt (1987). 



model and when n = 1 it is equivalent to an exponential 
model. The parameter 6„ has been approximated by b n — 
1.9992n-0.3271, for 0.5 < n < 10 (Capaccioli 1989; Prugniel 
& Simien 1997). The exponential model for the disc surface 
brightness profile can be written as follows, 



n(R) = no exp (—R/h), 



(3) 



where is the disc central surface brightness and h is the 
disc scale length. Our bulge-disc decomposition ignores the 
inner 4" of M31, which is dominated by an independent nu- 
clear feature (Light, Danielson & Schwarschild 1974). The 
results of the bulge-disc decomposition can be seen in Figure 
[1] Note that although there appears to be a slight break in 
the disc surface brightness profile at ~4000", the residuals 
are still very small. This break may appear due to the pres- 
ence of spiral structure, evident at similar radii in the 24 /jm 
of M31. Observational properties of M31 are derived from 
this bulge-disc decomposition and these are listed in Table 
1. 

With a disc scalelength of 25.92' and a surface bright- 
ness profile extending to a radius of 1?5, we are tracing the 
disc profile out to ~3.5 scalelengths. Assuming that the ex- 
ponential disc profile extrapolates to infinity, our profile in- 
cludes 86.1 per cent of the total disc light. Below, where we 
assign a mass to the disc, we extrapolate the disc profile to 
6 scalelengths, which would contain 98.2 per cent of the disc 
light. Since Ibata et al. (2005) have shown that the disc of 
M31 extends to very large radii (~70 kpc) with the same 
scalelength, our above calculation assumes no disc trunca- 
tion. Both the 10 kpc ring and the 2.6 kpc bar (Beaton et al. 
2007) show up as small perturbations in the surface bright- 
ness profile (Fig. 1) at radii of 46.67' and 26.67' respectively. 

We now assign masses to the disc and bulge of M31. 
The stellar mass-to-light ratio in the JG-band is a well cal- 
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Figure 3. B — R colour profile for M31 from Walterbos & Ken- 
nicutt (1987) corrected for foreground Galactic extinction. 



ibrated quantity (see Bell et al. 2003) which depends on 
B — R colour. However, a similar calibration does not ex- 
ist for the Spitzer 3.6 pm waveband. An alternative method 
would be to estimate the iG-band surface brightness profile, 
given the 3.6 /im profile, and since we know the B — R colour 
profile, we have chosen to determine the K— [3.6] colour de- 
pendence on B — R colour using both the Bruzual & Chariot 
(2003; hereafter BC03) stellar population synthesis models 
and the Maraston (2005; hereafter M05) stellar population 
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Table 1. M31 Observational Data. Hubble type taken from de Vaucouleurs et al. (1991). Distance in kpc taken from Holland et al. 
(1998). 



Parameter 



Measurement 



Hubble type 
Distance (kpc) 

Position Angle of major axis (degrees) 
Bulge effective radius, R e (arcmin) 
Bulge effective radius, R e (kpc) 

Bulge surface brightness at the effective radius, /t e (3.6 fim-mag arcsec - 2 ) 
Bulge Sersic index, n 

Disc central surface brightness, /to (3.6 /im-mag arcsec -2 ) 

Inclination corrected disc central surface brightness, fio i (3.6 /im-mag arcsec -2 ) 
Disc scalelength, h (arcmin) 
Disc scalelength, h (kpc) 
Disc luminosity, L disc (Lq) 
Bulge-to-disc ratio, B/D 



SA(s)b 

784 

45 

8.45±0.43 

1.93±0.10 

16.42±0.83 

1.71±0.11 

17.35±0.93 

18.94±0.93 

25.92±1.28 

5.91±0.27 

(6.08±0.60)xl0 1 

0.57±0.02 



synthesis models. The BC03 models use the stellar library of 
Pickles (1998), which includes stars in a wide range of spec- 
tral types (05-M10) and luminosity classes (I-V) and three 
metallicity groups (metal poor, metal rich and solar). The 
Pickles (1998) library has wavelength coverage from 0.1 )im 
to 2.5 /im. BaSeL 3.1 spectra (see BC03 Tables 2 and 3) are 
then used to extend the Pickles (1998) spectra into the in- 
frared from 2.5 /im to 160 /im. The spectra of M giant stars 
are the only ones not observed in the Pickles (1998) library. 
As a result the BC03 models use synthetic M giant spectra 
calculated by Fluks et al. (1994). These spectra are extended 
into the infrared using the model atmospheres of Schultheis 
et al. (1999), which cover a spectral range from 0.5 /im to 10 
/im, and 10 equally spaced stellar temperatures in the range 
2600 < T cff < 4400 K. The main difference between the 
BC03 and the M05 models is the way in which thermally- 
pulsing AGB (TP-AGB) stars are treated. The M05 models 
provide a better treatment of TP-AGB stars than the BC03 
models. Furthermore, Maraston et al. (2006) find that in 
general the M05 models provide better fits to observations 
of infrared spectra of galaxies than the BC03 models. They 
indicate systematically lower ages and, on average, masses 
that are 60 per cent lower for the stellar population samples 
in these galaxies. 

Figure [2] shows the result of both the BC03 (left panel) 
and M05 (right panel) models. In both cases, the dotted 
lines in Figure [5] represent the B — R colour of the inner 
and outer disc of M31 (Walterbos & Kennicutt 1987 and 
Figure [3| after correction for Galactic foreground extinc- 
tion. Since the M05 models provide an updated treatment 
of TP-AGB stars and generally provide better fits to ob- 
servations than the BC03 models, from here on we con- 
centrate on the M05 results. However, it is interesting to 
note that the if— [3.6] colour index predicted by the BC03 
and M05 models differs by <0.03 mag (see Figure 2). Lines 
of constant metallicity (solid lines) and age (dashed lines) 
are shown. The typical metallicity in the inner few kpc of 
M31 is ~2.4 Zq, whereas in the outer disc this falls to ~0.4 
Z Q (Brewer, Richer & Crabtree 1995; Worthey et al. 2005). 
From this metallicity gradient and the M05 models (Figure 
[5] right panel) it can be seen that the B — R colour gradi- 
ent in M31 is almost entirely due to changes in metallicity, 
and that the A"— [3.6] colour difference is < 0.1 mag. Fur- 



thermore, we can now use the B — R gradient (Figure [3| 
and the metallicity gradient of M31 to estimate the K— [3.6] 
colour as a function of radius. We find that the inner re- 
gion of M31 has a colour K— [3.6] ~ 0.25, and the outer disc 
(at 26 kpc) has a colour K-[3.6]~ 0.18. Given the B - R 
colour profile and the K-band mass-to- light ratio calibration 
from Bell et al. (2003) we calculate a central mass-to-light 
ratio, M/Lk = 1.0 ± 0.1 and a mass-to-light ratio at 26 
kpc, M/L K = 0.75 ± 0.1. From the K-[3.6] colours de- 
termined in our stellar population synthesis analysis, the 
above mass-to-light ratios are equivalent to 3.6 /im mass- 
to-light ratios of M/L3.6 — 1.25 ± 0.10 in the centre and 
M/L3.6 ~ 0.89 ± 0.10 at 26 kpc (equivalent to a M/L gra- 
dient 0.13 kpc - ). In our analysis, we prefer to keep the 
central (M/Lz.§) value a free parameter, and as such we al- 
low it to vary in the range 0.05 < (M/L3.6) < 1.5. This is 
a standard approach to adopt in mass modelling of galax- 
ies (see Humphrey et al. 2006; Zappacosta et al. 2006) and 
as discussed at the end of section 3.2 and shown in Table 
2, when the M/L3.6 is left as a free parameter, the best- 
fit value found is M/L3.6 = 1.15, i.e. consistent with the 
value of M/L3.6 derived here. We use our derived value for 
the 3.6 /im disc luminosity of L3.6 = (6.08 ± 0.60) x 10 10 
Lq. We then use the derived disc and bulge light profiles 
L3.6 = I/disc + ibuige to determine the stellar mass contribu- 
tion to the rotation curve, M* = (M / X3. 6)^3.6- 

A concern in using the 3.6 /im Spitzer waveband to de- 
termine the underlying stellar mass, is the effect of emission 
from hot dust in this waveband, although this is probably 
only important in or near HII regions. In order to place some 
constraint on this, we have chosen to explore the emission 
from dust in the the near-infared A'-band at 2.2 /im. Us- 
ing near-infrared spectroscopy at 2.2 /im it has been shown 
that hot dust can account for up to 30 per cent of the contin- 
uum light observed at this wavelength in areas of active star 
formation, i.e., spiral arms (James & Seigar 1999). When 
averaged over the entire disc of a galaxy, this reduces to a 
<2 per cent effect, if one assumes that spiral arms can be 
up to 12° in width. At 3.6 /im this would therefore result in 
<3 per cent of emitted light from dust. 

Another concern for the 3.6 /im waveband, would be 
the contribution from the Polycyclic Aromatic Hydrocarbon 
(PAH) emission feature at 3.3 /im. However, ISO spectra 
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covering the wavelength range 5.15 to 16.5 /im have shown 
that PAH emission in M31 seems to be very weak (e.g., Ce- 
sarsky et al. 1998). Furthermore, in a survey of actively star 
forming galaxies, Helou et al. (2000) found that the 3.3 /im 
feature was also very weak when they analysed the average 
2.5 — 11.6 /im spectrum of 45 galaxies. The contribution of 
the PAH feature to the 3.6 /im Spitzer waveband, is there- 
fore not a major concern. 

In the following section, where we model the dark mat- 
ter halo, we have included a bulge surface brightness model 
using a deprojected Sersic function. The function we use is 
from Trujillo et al. (2002) who show that the deprojected 
mass density profile is given by, 

P{R) = R~^ l=C(h) (4) 

where R e is the effective radius, 1(0) is the central in- 
tensity, T is the stellar mass-to- light ratio, h = R/R e , 
C(h) = hi(\ogh) 2 + h 2 \ogh + h 3 and K v is the i/ th -order 
modified Bessel function of the third kind. Trujillo et al. 
(2002) give values for u, p, hi, /12 and /13 for different values 
of the Sersic parameter, n. In our case, where n — 1.71, we 
adopt values of v = 0.45394, p = 0.67930, hi = -0.0162935, 
hi = -0.18161 and h 3 = 0.04435. Finally, the constant f 1/2 
depends on the three-dimensional spatial orientation of the 
bulge. In our case, we assume the bulge is intrinsically spher- 
ical, and therefore f 1 ^ 2 = 1. 

It should be noted that this version of the deprojected 
Sersic function (equation [4} is an approximation, but it is 
good enough for our purposes. 

3.2 Modelling the dark matter halo 

A range of allowed dark matter halo masses and density 
profiles are now explored, using two models for disc galaxy 
formation. In the first we assume that the dark matter halo 
surrounding M31 has not responded significantly to the for- 
mation of the disc, i.e., adiabatic contraction does not occur 
(the "no- AC" model). In this case the dark matter contri- 
bution to the rotation curve is described by a density profile 
similar to those found in dissipationless dark matter simu- 
lations (the NFW profile): 

P{R) = (R/R s )(l + R/R s y (5) 

where R s is a characteristic "inner" radius, and p s is a cor- 
responding inner density. This is a two parameter function 
and is completely specified by choosing two independent pa- 
rameters, e.g., the virial mass M vlv (or virial radius R V1I ) and 
concentration c v i r = R v i T /R 3 (see Bullock et al. 2001a for a 
discussion). Similarly, given a virial mass M v j r and the dark 
matter circular velocity at any radius, the halo concentra- 
tion Cvir is completely determined. 

In the second class of models we adopt the scenario of 
adiabatic contraction (AC) and specifically use the origi- 
nal scheme discussed by Blumenthal et al. (1986, hereafter 
B86; see also Bullock et al. 2001b and Pizagno et al. 2005). 
We also investigate the slightly revised scheme presented in 
Gnedin et al. (2004, hereafter G04). As we discuss below, 
the G04 scheme reproduces the outer slope of the observed 
rotation curve of M31, but needs a very high concentration 
to do so. 



For each AC algorithm (B86 and G04) and for the no- 
AC model we generate a grid of final rotation curves. We 
vary the baryonic contribution to the rotation curve by al- 
lowing the bulge-disc decomposition parameters for h and 
Ldisc to range over their ±2<r values (see Table 1 and the pre- 
vious section) . The central mass-to-light ratio of the galaxy 
is allowed to vary over the range 0.05 < (M/L) < 1.5. From 
this we determine the bulge and disc contribution to the 
rotation curve. For the bulge a simplified spherical model 
is used. For the disc the potential is derived assuming zero 
thickness and modified Bessel functions as described by Bin- 
ney & Tremaine (1987). For each set of baryonic parameters 
we generate a range of dark halo models with total circular 
velocities at 2.2 disc scalelength (V2.2) that span the best-fit 
value for M31 (V2.2 = 270 km s" 1 ) quoted by Rubin & Ford 
(1970) by a wide margin: 200 < V2.2 < 340 km s" 1 . In prac- 
tice we achieve each V2.2 value in our models by allowing the 
initial halo concentration to vary over > 4cr of its expected 
range c v i r = 3 — 60 (Bullock et al. 2001a; Wechsler et al. 
2002; Maccio et al. 2006) and then setting the halo virial 
mass M v i r necessary to reproduce the desired value of V2.2- 
Finally, we determine the implied fraction of the mass in 
the system in the form of stars compared to that expected 
from the Universal baryon fraction, /* = M*/(/(,M v i r ) and 
demand that /* obeys 0.01/b < /* < /{,. From this grid of 
choices we derive best-fitting dark halo parameters by min- 
imizing the reduced-x 2 for our two choices of rotation curve 
data. These results are summarized in Table 2. 

For the most part, when minimizing \ 2 , we assume that 
the errors on the observational data (i.e., the rotation veloc- 
ities) are not correlated. Instead we assume that they are 
Gaussian in nature. However, this is not the case for the 
Carignan et al. (2006) HI data (see Figure U top row and 
bottom right beyond 25 kpc), where the errors obviously are 
correlated. In order to account for the strange behaviour of 
the errors on these HI rotation velocities we have treated 
all of this data as an average rotation velocity of 226.3 km 
s _1 at 31.16 kpc. This is an idealized simplification but al- 
lows a quantitative estimate of the relative goodness of fit 
between models. Furthermore, although error bars are not 
available for the "M2" curve, we assume 10 per cent errors 
in the model fitting, as this is similar to the size of the er- 
rors quoted for the Rubin & Ford (1970) Ha rotation curve. 
We therefore quote a x 2 f° r the "M2" curve in Table 2 that 
is consistent with 10 per cent errors. We then determine 
how consistent all the derived rotation curves are with each 
other, by adopting the formal best fitting rotation curve as 
our fiducial mass model. For our "Ml" curve, the number 
of degrees of freedom in the fit is v = 31. The number of 
data points we fit to is n^ata = 37, and 6 parameters are de- 
termined from the fit, the NFW concentration c v j r , the halo 
virial mass M v i r , the stellar mass-to-light ratio M/L, the 
disc luminosity Ldisc, the disc scalelength h and the bulge- 
to-disc ratio B/D, thus leaving the degrees of freedom in 
the fit as v — 31. 

The Ml rotation curve from Rubin & Ford (1970) is 
shown by the solid squares with error bars in Figure |4]( top 
left, top right, middle right, bottom left) along with our best- 
fitting rotation curve models (solid lines). The contribution 
from stellar mass (short-dashed) and dark matter (long- 
dash) in each case is shown along with the rather minor con- 
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Figure 4. Top left: Rotation curve data (solid square points) from Rubin & Ford (1970) with the best fit model, with adiabatic contraction 
(Ml B86 model), overlaid (solid line). For comparison the Hi rotation velocities (grey open triangles) from Carignan et al. (2006) are 
shown (the black point corresponds to the average velocity and radius of the Carignan data that we make use of). The best fit rotation 
curve is decomposed into four components, the contribution from dark matter (short-dashed line), the contribution from baryonic matter 
(long-dashed line), the Hi contribution (dotted line) from Carignan et al. (2006) and the contribution from the central 10 s Mq black 
hole (dot-dashed line) from Bender et al. (2005). Middle left: Rotation curve CO and Hi data (solid square points) from Loinard et al. 
(1995) and Brinks & Burton (1984) respectively, as used by Klypin et al. (2002). The best fit model with adiabatic contraction (M2 B86 
model) is overlaid (solid line). It is decomposed into four components as in the top left panel. Bottom left: Rotation curve data from 
Rubin & Ford (1970; solid points) and Carignan et al. (2006; triangles) and the data as used by Klypin et al. (2002; hollow circles), with 
adiabatic contraction model (M1+M2 B86), overlaid. Top right: The same as the top left panel, but with the best fitting model using a 
modified version of adiabatic contraction (Ml G04 model) as in Gnedin et al. (2004). The long-dashed curve includes the total baryonic 
mass, including the Hi mass and the central black hole. Middle right: The same as the top left panel, but with the best fitting model 
without an adiabatically contracted halo (Ml no-AC model). In all five cases the NFW concentration, c v i r , and virial mass, M v - lr , are 
shown. 
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Table 2. M31 best fitting models. "Ml B86" is the best-fit model to the Ha rotation curve from Rubin & Ford (1970) (Ml) using the 
AC prescription of B86. "M2 B86" is the best-fit B86 model to a combination of the CO rotation velocities from Loinard et al. (1995) 
and the Hi velocities from Brinks & Burton (1984) (M2). "M1+M2 B86" is the best-fit B86 model to the Ha rotation curve from Rubin 
& Ford (1970) (Ml) combined with the CO rotation velocities from Loinard et al. (1995) and the Hi velocities from Brinks & Burton 
(1984) (M2). "Ml G04" corresponds to the Ml data fit using the modified adiabatic contraction model of Gnedin et al. (2004) and "Ml 
no- AC" , uses the same data without adiabatic contraction. 



Parameter 



Ml 

B86 



M2 
B86 



M1+M2 
B86 



Ml 
G04 



Ml 
no-AC 



Shear 

V2.2 (km s" 1 ) 

Disc luminosity, Ldisc (Lq) 

Disc scale length, h (kpc) 

Central disc mass-to-light ratio, 
M/L (3.6 /im solar units) 

Mass-to-light ratio gradient, 
d{M/L)/dR (kpc" 1 ) 

Initial NFW concentration, c v i r 

Bulge-to-disc ratio, B/D 

Virial mass, M v - lr (Mq) 

Dark matter concentration, 
C DM> R< 10 kpc 

Total mass concentration, 
ctoti R < 10 kpc 

Stellar baryon fraction, /* 
X 2 

Degrees of freedom, v 



0.51 0.50 0.51 0.54 0.50 

275±5 275±5 275±5 275±5 275±5 

(4.9 ±0.6) x 10 10 (6.1 ± 0.6) x 10 10 (5.5 ± 0.6) x 10 10 (4.9 ± 0.6) x 10 10 (6.1 ± 0.6) x 10 10 

5.09±0.06 5.91±0.06 5.50±0.06 5.50±0.06 5.91±0.06 

1.15±0.10 1.00±0.10 1.05±0.10 1.05±0.10 0.95±0.10 



0.013 



0.013 



0.013 



0.013 



0.013 



20.0±1.1 16.0±1.1 18.0±1.1 40.0±1.1 51.0±1.1 

0.57±0.02 0.57±0.02 0.57±0.02 0.57±0.02 0.57±0.02 

;.2 ± 0.2) x 10 11 (8.5 ± 0.2) x 10 11 (8.4 ± 0.2) x 10 11 (8.9 ± 0.2) x 10 11 (7.3 ± 0.2) x 10 11 
14.9±1.4% 14.4±1.5% 14.7±1.5% 9.9±1.3% 14.1±2.0% 



26.2±2.5% 

56.7±3.9% 
38.75 
31 



25.4±2.4% 

59.2±4.5% 
40.63 
27 



25.9±2.4% 

58.1±4.2% 
80.31 
64 



23.8±2.7% 

56.3±6.0% 
46.50 
31 



26.5±2.6% 

63.5±5.6% 
41.58 
31 



tributions from neutral gas (dotted) and the central black 
hole (dot-dashed). 

Overlaid in the top left panel of Figure [4] is a model that 
includes adiabatic contraction (overall reduced-^ 2 = 1.25). 
The dark halo in this case has undergone adiabatic con- 
traction according to the B86 prescription. Our best-fitting 
parameters are summarized in Table 2 under "Ml B86". We 
determine the halo of M31 to have a total (dark+baryonic) 
virial mass of M vir = (8.2 ± 0.2) x 1O U M and an initial 
concentration c v ir = 20.0 ± 1.1. Our derived virial mass is 
almost a factor of 2 smaller than that quoted by Klypin et 
al. (2002). We also find a concentration that is higher than 
both Klypin et al. (2002), who find c v i r = 12 and Widrow et 
al. (2003), who find c v i r = 10. In the case of Widrow et al. 
(2003) it should be noted that they initially fit an isothermal 
model given by Evans (1993). They then fit an NFW profile 
to this model. Given that the isothermal model favors a con- 
stant density core, the derived NFW concentration will be 
articifially low and fitting an NFW model in a more direct 
manner may result in a higher concentration. 

The model rotation curve fits the data adequately 
within the uncertainties. The innermost data point from the 



observational rotation curve is significantly lower than the 
model. This can be explained by the recent strong evidence 
for a bar, which dominates the near-infrared light out to a 
radius of ~2.6 kpc, in M31 (Athanassoula & Beaton 2006; 
Beaton et al. 2007). This bar has been used to explain the 
minimum seen in the observed rotation curve at a radius of 
~ 2.7 kpc (Athanassoula & Beaton 2006). This bar may also 
affect the outer rotation velocities as well, although proba- 
bly by less than 10 per cent. This would therefore be within 
the errors associated with the observed data. Indeed, shift- 
ing the outer rotation velocities systematically upward by 
~10 per cent changes the best-fit models by an amount that 
remains within the quoted error range of the original fit (see 
Table 2). It should also be noted that by a radius of ~ 26 — 35 
kpc the rotational velocity has dropped to ~ 225 — 230 km 
S _1 , which is similar to the rotational velocities seen at radii 
of 25 — 35 kpc in the HI rotation curve (Carignan et al. 2006; 
open triangles). 

The best-fitting model without adiabatic contraction is 
shown by the solid line in the middle right panel of Figure 
[4] This, the no-AC model, was determined in exactly the 
same way as our best-fit B86 AC model. As listed in Ta- 
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Figure 5. Virial mass as a function of NFW concentration from 
our modelling with our la, 2a and 3a confidence contours shown 
as solid lines for our "Ml B86" model and as dashed lines for our 
"Ml no- AC". 



ble 2, the no- AC assumption requires a much higher NFW 
concentration, c v i r = 51.0 ± 1.1 and a lower virial mass, 
Mvir = (7.3 ± 0.2) x 1O U M . This model also produces 
a reasonable fit, with a reduced-;^ 2 = 1.35, but it should 
be noted that the required concentration is, at best, only 
marginally acceptable within a ACDM context. This would 
correspond to a ~ 4.4er outlier from the expected concen- 
tration distribution of M31-sized haloes (e.g., Bullock et al. 
2001a). 

In Figure [5] we plot the NFW concentration, c v i r , as a 
function of the virial mass, M V ir, for both the "Ml B86" 
case (solid line) and the "Ml no-AC" case (dashed line). 
The contours represent the la, 2a and 3cr confidence lev- 
els for the concentration and virial mass for both models. 
From figure [5j one can see that the "Ml no-AC" model is 
inconsistent with any reasonable concentration of c v ir < 30 
at the 5ct level at least. However, since it is possible to pro- 
duce a reasonable fit with a pure NFW model, albeit with 
a very large concentration, we do not rule out our "Ml no- 
AC" case, although it is clearly not consistent with standard 
ACDM cosmology. 

The confidence contours shown in Figure [5] have also 
allowed us to estimate errors on the derived best-fit masses 
and concentrations. The la, 2a and 3a contours are drawn 
as the boundaries where x 2 increases by 1.0, 2.71 and 6.63 
relative to its minimum value. For our best-fit 'Ml B86' 
model, x 2 — 38.75. There is also a 'Ml B86' model for which 
the x 2 =39.75, i.e. ler higher than the best fit. For this model, 
the concentration is c v ir=21.1 and the virial mass is M v i r = 
8.0 x 10 Mq. There is a change in the concentration of 
Ac v ir=l.l and a change in the virial mass of AA/ V i r = 0.2 x 
10 11 M© at the la level, i.e. when the value of \ 2 is increased 
by 1. The values of AcVir and AM v j r are therefore adopted 
as our la error on c v i r and M v j r respectively, and these are 
the errors that we report in Table 2. 



An intermediate result is shown in the top right panel 
of Figure [4] Here we use the AC model of G04, which pro- 
duces less contraction than the B86 prescription. If the ini- 
tial NFW concentration is held fixed at the same value as 
the best-fitting B86 case, i.e., c v j r = 20, then the G04 model 
overestimates the observed rotation curve at r > 15 kpc. 
The G04 model can only reproduce this outer slope when 
the initial halo has c v j r = 40. With this high concentration 
the reduced-x 2 of the G04 model is 1.50. This model can- 
not be ruled out, but (as in the pure NFW case) its halo 
concentration is higher than expected in ACDM (Bullock et 
al. 2001a, Maccio et al. 2006). Therefore, from here on, we 
chose to focus on the B86 AC case, but provide the best- 
fitting model results for the G04 case in Table 2. 

For completeness, we also perform a fit to the same ob- 
served rotation curve that was used by Klypin et al. (2002) in 
their M31 mass decomposition. Our best fitting mass model 
is shown in the middle left panel of Figure [4] ( "M2 B86" in 
Table 2). Here we have adopted the B86 AC model. The 
initial NFW concentration in this case (c v i r = 18) is higher 
than that the value of c v i r = 12 favored by Klypin et al. 
(2002) and lower than our "Ml B86" case. However, it is 
not a large difference, and at face value our "Ml B86" and 
"M2 B86" models are consistent with each other. As we 
discuss below, the difference between our result and that 
of Klypin et al. (2002) is driven mainly by our updated 
baryonic model, which seems to provide a less concentrated 
bulge. For the M2 case, we calculate a x 2 , assuming 10 per 
cent errors. As a final consistency check, we also perform a 
fit of the B86 AC model to a combination of the Ml and M2 
data (bottom left of Figure g] "M1+M2 B86" in Table 2). 
This produces a result consistent with "Ml B86" and "M2 
B86". 

Figure |S] shows the enclosed mass as a function of ra- 
dius for our fiducial model (Ml B86). The fraction of mass 
contained within the central 10 kpc (the "central mass con- 
centration" hereafter, as originally defined by Seigar et al. 
2006) is ctot = 26.2 ± 2.5 per cent. In absolute terms, the 
mass contained within 10 kpc is (21.5 ±2.6) x 10 10 M Q . This 
is about a factor of 2 larger than that calculated by Rubin 
& Ford (1970), but it should be noted that they did not 
take into account non-circular orbits due to the presence of 
a bar, which we now know exists in M31. This would have 
the effect of increasing the mass contained within the bar 
region. The mass calculated within a 35 kpc radius from the 
HI rotation curve by Carignan et al. (2006) is 3.4 x 10 11 Mq, 
similar to the mass derived from our best fit model within 
a 35 kpc radius of (3.4 ± 0.5) x lO n M . Within a 31 kpc 
radius we find a mass of (3.1 ± 0.4) x W 1x Mq, which is sim- 
ilar to the mass of 2.8 x IO^Mq within the same radius 
found using kinematic data of planetary nebulae (Evans & 
Wilkinson 2000). The total implied stellar baryon fraction 
of this model is /* = 56.7 ± 3.9 per cent. 

Other mass models of M31 (Klypin et al. 2002; Geehan 
et al. 2006; Widrow et al. 2003) have used bulge velocity dis- 
persions as inputs into their models. For example, Klypin et 
al. (2002) used the line-of-sight velocity dispersion of M31 
from Kormendy & Bender (1999) to constrain mass mod- 
els of M31 within the bulge region. In the bulge region, the 
velocity dispersion, a is dominant (c. 150 km s _1 ) at least 
out to twice the effective radius, i.e. 2R e =3.96 kpc. Gee- 
han et al. (2006) assume a Hernquist (1990) model for the 
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Figure 6. Enclosed mass as a function of radius for model Ml 
B86. The solid line indicates the total mass. Also indicated are 
the HI mass out to a 35 kpc radius (dotted line), the bulge+BH 
mass (short-dashed line), the disc mass (long-dashed line) and 
the halo mass (dot-dashed line). The data points indicate masses 
derived using other observational methods at fixed radii of 10 kpc, 
31 kpc, 35 kpc and 125 kpc from Rubin & Ford (1970), Evans & 
Wilkinson (2000), Carignan et al. (2006) and Fardal et al. (2006) 
respectively. 

bulge, and compare the expected velocity dispersion pro- 
file from Jeans equation modelling with the observations of 
Kormendy (1988). Our approach differs from this. We have 
chosen to assume a Sersic law for the Bulge and we then 
perform a mass-to-light conversion. In our approach, we es- 
sentially ignore any kinematics within the inner few kpc, 
because in that region the kinematics are likely to be per- 
turbed by the bar or affected by streaming motions along the 
bar. Instead, we assume that stellar light traces stellar mass. 
Figure 4 shows that the bar kinematics are strongest within 
the inner ~6 kpc. Beyond 10 kpc, the bar kinematics are not 
likely to affect the disc kinematics. Furthermore, we assume 
that the bulge is spherically symmetric (an assumption also 
made in the papers by Klypin, Geehan and Widrow). The 
bulge in our model has a mass of Mbuige — 3.5 x 10 10 Mq. It 
is reassuring to note that the Jeans modelling of Geehan et 
al. (2006) derives a bulge mass of M bu i ge = 3.2 x 10 10 M 0) 
very similar to our bulge mass. 

We note that the bulge component in our models is 
less concentrated than that derived by some other groups. 
Klypin et al. (2002) have a more concentrated bulge, but 
that is probably due to a combination of the fact that they 
use a Hernquist model (equivalent to a Sersic index of n — 4) 
for the M31 bulge and optical data to determine the bulge 
characteristics. Geehan et al. (2006) also use a Hernquist 
model. This type of model will naturally produce a more 
concentrated bulge than a Sersic model with index n ~ 1.7. 
The use of the M/L ratios from Bell & de Jong (2001) also 
gives Klypin et al. (2002) an overall more massive bulge. 
The more recent models by Tamm et al. (2007) and Tempel 
et al. (2007) make use of a model that is similar to a Sersic 



law, but they find an index of n ~ 3 compared to our index 
of n ~ 1.7. As a result, their bulge model is also naturally 
more concentrated than ours. It is also interesting to note 
that replacing our Sersic bulge with a bulge profile modelled 
using a Hernquist or exponential law does not affect the 
overall dark matter halo profile by more than 10 per cent. 
This may be because the overall bulge mass does not change; 
It is just the way in which the bulge mass is distrbuted that 
changes. 

The issue of turbulent gas motions is very important for 
low-mass galaxies with small rotation velocities (e.g., Valen- 
zuela et al. 2007), as this requires a correction for pressure 
support or asymmetric drift in the mass model. The cor- 
rection to the total mass is ~ (e>"gas/Kot) 2 . Since the Ha 
and HI rotation curves are measured from gas, it is impor- 
tant to understand what effect the turbulent gas motions 
may have on our mass model. Since, M31 is a galaxy with 
a high rotation velocity (~ 275 km s _1 ) and a turbulent 
velocity dispersion in galaxy discs is typically < 30 km s _1 
(i.e. around 10 per cent of the rotation velocity) this effect is 
small. Indeed, if we assume that the scatter of the Ha Rubin 
& Ford (1970) rotation curve around our best fitting model 
is due to turbulent gas motions, then we can derive an esti- 
mate for their typical value in M31. We can do this because 
the Rubin & Ford (1970) rotation curve is derived from Ha 
spectroscopic observations of individual HII regions in M31. 
In this case, we have to assume that the errors on the Ru- 
bin & Ford (1970) rotation curve are uncorrelated, and we 
derive a turbulent velocity dispersion of (J gas = 18.3 km s _1 . 
Adding this in quadrature to the maximum rotation velocity 
from the Rubin & Ford (1970), results in a mass difference 
of ~0.4 per cent for the total (baryonic + halo) galaxy mass. 

As mentioned above, the preferred total 
(dark+baryonic) virial mass for M31 is (8.2 ±0.2) x 10 n M Q . 
This is similar to the estimate of ~ (7.9 ± 0.5) x lO n M 
derived from satellite kinematics (Cote et al. 2000) and 
the total mass of ~ 8 x lO n M0 from kinematics of the 
Andromeda Stream (Fardal et al. 2006; Geehan et al. 2006). 
It is also close to the lower limit of 9 x 10 Mq derived from 
kinematics of halo stars (Chapman et al. 2006) . Our favored 
virial mass is about a factor of 2 smaller than the value 
favored by the Klypin et al. (2002) analysis of the M31 
rotation curve. This can be attributed to three differences in 
our analysis. First, Klypin et al. (2002) use the HI rotation 
curve of Brinks & Burton (1984) to find their best fitting 
model. In this study we use the Ha rotation velocities of 
Rubin & Ford (1970). However, we find that this is not the 
major driver. We find that the model rotation velocities 
are consistent with those from the more recent HI study 
of Carignan et al. (2006), and when we fit to the same 
observed rotation curve as Klypin et al. (2002), we find that 
our best-fit model is consistent with our best-fit model to 
the Rubin & Ford (1970) curve. Furthermore, if we try and 
model a combination of data from Rubin & Ford (1970) 
and Klypin et al. (2002) we also find a best-fitting rotation 
curve consistent with our fiducial case. More importantly, 
Klypin et al. (2002) used the M/L ratios of Bell & de Jong 
(2002) to determine the stellar mass from their bulge-disc 
decomposition. In this study we allow M/L to vary freely 
in the range 0.05 < (M/L) < 1.50, and find a best-fitting 
value close to the expected 3.6 /im mass-to-light ratios (see 
section [3~T) . Possibly the most significant difference is that 
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Klypin et al. (2002) perform a bulge-disc decomposition 
using an R band image, whereas we use a Spitzer 3.6 /im 
image. Using a 3.6 ^m image provides a better trace of 
the stellar mass (e.g., Barmby et al. 2006 and references 
therein). Overall, we should therefore have a more accurate 
estimate of the baryonic mass profile in M31. Indeed, if we 
adopt the baryonic mass profile of Klypin et al. (2002) and 
try to fit either to the HI or to the Ha rotation curve, we 
find a best fitting model that is consistent with their result. 

Typically, the best-fitting M/L values we find are in 
the range 0.95-1.15, depending on the model used. The low- 
est value of M/L — 0.95 is difficult to reconcile with the 
expected value, which is M/L3.6 = 1.2 ± 0.1, for a galaxy 
with the central B — R colour of M31. This expected value 
was determined by combining the Bell et al. (2003) A^-band 
M/L with the K— [3.6] colours derived from our stellar pop- 
ulation synthesis described in section 3.1 (Figure [2}. The 
lowest value we find (M/L = 0.95) is for our no-AC model. 
We also find that this model is difficult to explain in the con- 
text of ACDM cosmology, due to its large concentration. Our 
Ml B86 model reproduced the highest value of M/L = 1.15 
that we find. This is consistent with the expected range of 
M/L values. 



4 SUMMARY 

In this paper we have derived new mass models for M31, 
based on an updated baryonic contribution to the rotation 
curve. The baryonic contribution is derived using a Spitzer 
3.6 /im image. We find that both a pure NFW and an 
adiabatically contracted NFW profile can produce reason- 
able fits to the observed rotation curve of M31. However, 
the best-fit concentration of M31 for a pure NFW model is 
Cvir = 51 ± 1.1, which is disfavoured when compared to the 
expected value of log 10 (c v ir) = 1-08 ± 0.14 (Wechsler et al. 
2002; Maccio et al. 2006) in a ACDM cosmology (our best 
fitting concentration is a 4.4a outlier). The best fit produced 
by the AC recipe of Gnedin et al. (2004) brings down the 
concentration to c v i r = 40, or log 10 (c v i r ) = 1.60, and is thus 
marginally consistent with ACDM as a ~3.7<r outlier. If one 
accepts ACDM cosmology, then the most consistent fit to 
the observed rotation curve seems to be that produced by 
the AC recipe of Blumenthal et al. (1986). In this case the 
concentration, c v i r = 20 ±1.1, and this is consistent with the 
expectations of ACDM, to within 1.5a. 

For all of the types of model we use, the derived best-fits 
are in good agreement with the mass distributions derived 
from various other observational methods. Our estimate of 
the halo virial mass lies in the range M v i r = (7.3 ± 0.2) x 
lO n M (for the Ml no-AC case) and M vil = (8.9 ± 0.2) x 
10 n M© (for the Ml G04 case). This is in close agreement 
with the virial mass of ~ 8 x 10 11 A/q found via kinematics 
of satellite galaxies (Cote et al. 2000) and the Andromeda 
Stream (Fardal et al. 2006; Geehan et al. 2006), and the 
lower limit of 9 x 10 11 Mq found from kinematics of halo 
stars (Chapman et al. 2006). Our Ml B86 model, which is 
in closest agreement with the expectations of ACDM, lies 
in the middle of the range with a virial mass of M v i r = 
8.2 ± 0.2 x 10 n M Q . 

Another recent M31 mass model by Tamm, Tempel & 
Tenjes (2007) and Tempel, Tamm & Tenjes (2007) also uses 



Spitzer 3.6 fim imaging to determine the total mass of M31, 
and they find a mass of l.lxl0 12 M Q for a pure NFW model, 
i.e. a slightly higher mass than for our best fit. However, they 
did not look at the question of AC. Also, their bulge seems 
to have a more concentrated density profile (see Figure 6 
in Tempel et al. 2007), which is probably driven by the use 
of rotation velocities inside 5 kpc. This may be driving the 
differences between our model and the model presented in 
Tempel et al. (2007). We have been cautious not to include 
rotation data within the inner few kpc in our model, since it 
is unlcear how non-circular motions induced in this barred 
region will affect the overall mass model. 

It is interesting that, while it has been shown that adi- 
abatically contracted halo models do not generally fit the 
observed rotation curves of galaxies (e.g. Dutton et al. 2005, 
2007; Kassin et al. 2006a, b; Pizagno et al. 2005), the ob- 
served rotation curve of a well-studied galaxy such as M31 
can be well described by a halo that has undergone adiabatic 
contraction. Indeed, one either has to adopt (i) a pure NFW 
profile with an unphysically large concentration, within the 
context of ACDM, or (ii) an adiabatically contracted NFW 
profile and a more reasonable concentration. Such a scenario 
is in agreement with the results of Klypin et al. (2002), who 
claim that an adiabatically contracted halo is favoured for 
M31. 

Much of the power in this analysis came from including 
the normalization and outer slope of M31's rotation curve in 
our fits (where many similar studies of AC attempt to fit to 
just the rotation velocity at 2.2 disc scalelengths, V2.2). This 
highlights the usefulness of extended rotation curve data in 
constraining general models of galaxy formation. Having ro- 
tation curve data out to large radii is essential for determin- 
ing dark matter density profiles. 

Studies that test the applicability of adiabatic contrac- 
tion often make use of large samples of galaxies that are 
dominated by late- type disc-dominated galaxies, with small 
bulges or no bulge at all (e.g. Pizagno et al. 2005; Dutton 
et al. 2005, 2007; Gnedin et al. 2006; Kassin et al. 2006a, 
b). These "bulgeless" galaxies also tend to have very flat 
(or even rising) rotation curves, and it is therefore easy to 
explain these galaxy rotation curves without any need for 
adiabatic contraction. However, we suggest that early-type 
disc galaxies (e.g. Hubble types Sa-Sb), which have bulge- 
dominated centres, and rotation curves that fall-off due to 
the presence of a large bulge, may require adiabatic contrac- 
tion. Although we have studied very few galaxies so far (M31 
in this paper and two other galaxies in Seigar et al. 2006) 
our results indicate that in cases where rotation curves have 
a fall-off (two out of three cases), adiabatic contraction is 
favoured. This may be true for all galaxies with large bulges. 
If so, it may favour a secular (rather than merger-driven) 
origin for these bulge-dominated systems, as the gradual ac- 
cumulation of central mass increases the likelihood that AC 
will operate. 

It should also be noted that in our best-fit model, while 
it is consistent with the current HI data from Carignan et 
al. (2006), it is unclear if the rotation curve is declining or 
remaining flat past 30 kpc, due to the large error bars. The 
most important requirement for future kinematical studies 
of galaxies would be improved HI measurements at large 
radii. This will enable us to put better constraints on the 
halo properties of M31 and other spiral galaxies. 
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APPENDIX A: THE SPIRAL ARM PATTERN 
OF M31 

Because M31 is a highly inclined galaxy, it is difficult to learn 
how tightly wound its spiral arm pattern is from imaging 
data. However, Seigar, Block & Puerari (2004) and Seigar 
ct al. (2005, 2006) have shown that there is a very tight cor- 
relation between the spiral arm pitch angle (i.e. how tightly 
wound the spiral arms are) and rotation curve shear (mea- 
sured at a 10 kpc radius). Therefore the spiral arm pitch 
angle can be estimated from the shear measured for the M31 
rotation curve. Rotation curve shear is defined as 

c A 1 / R dV\ .. 

where A is the first Oort constant, u> is the angular velocity 
and V is the rotational velocity at radius R. The shear rate 
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depends upon the shape of the rotation curve. For a rotation 
curve that remains flat the shear rate, S = 0.5, for a falling 
rotation curve the shear rate, S > 0.5 and for a continually 
rising rotation curve the shear rate, S < 0.5. Following the 
method described by Seigar et al. (2006) we measure the 
shear at a 10 kpc radius. For M31, 10 kpc is equivalent to 
43.9'. The shear is measured using the prescription described 
by Seigar (2005) and Seigar et al. (2004, 2005, 2006). In 
this method an average slope is fit to the outer part of the 
rotation curve (i.e., past the solid-body rotation regime), 
and the shear is then calculated at the given radius. For M31, 
at a radius of 10 kpc, the shear is S — A/uj = 0.54 ± 0.02. 

There is little scatter in the correlation between shear 
rate and pitch angle, and the line of best fit is given by 

P = (64.24 ± 2.87) - (73.24 ± 5.53)5, (A2) 

where P is the spiral arm pitch angle in degrees and S is the 
shear. From this the pitch angle for M31 is P = 24?7 ± 4? 4. 
(It should be noted that the shear at 10 kpc determined for 
the Ml B86 model in section 2.2 is S = 0.52 and for Ml G05 
is S = 0.51. These estimates are consistent with the shear 
measured here. They result in pitch angles of P = 26? 2 ±2? 2 
and P — 26?8±2?2 respectively, both of which are consistent 
with the pitch angle derived above). 

Assuming a logarithmic spiral (see e.g., Seigar & James 
1998b for the equations that define a logarithmic spiral), 
the spiral arm pattern of M31 is calculated. Figure [T] shows 
the predicted spiral arm pattern of M31 for a face-on view 
(left) and with an inclination angle, i — 77? 5, and a position 
angle 45° (de Vaucouleurs et al. 1991) appropriate for M31 
(centre). It should be noted that the position and inclina- 
tion angles are not fitted parameters, they are taken from 
the literature. The spiral arm pitch angle calculated here is 
larger than the HI spiral arm pitch angle of 16° calculated 
by Braun (1991), although it is possible that the spiral arm 
pattern in the neutral gas may be a few degrees tighter than 
that observed in the underlying stellar distribution. A pitch 
angle of 16° would be a ~2cr outlier in the spiral arm pitch 
angle versus rotation curve shear correlation when the shear, 
S = 0.54 (see Figure 3 of Seigar et al. 2006). If this point 
was added to the correlation, it would be one of the larger 
outliers, yet the correlation would still be strong. For com- 



parison the 24 fim image is shown in the right panel of Figure 

m 

It is important to note that the spiral arm pattern 
shown in Figure [1] does not take into account the disc scale 
height and how broad the arms are. These plots are intended 
to illustrate show how tightly wound the arm pattern would 
be if viewed from face-on. When viewed close to edge-on (as 
is the case for M31), the scale height of the disc and breadth 
of the spiral arms make it virtually impossible to determine 
the spiral arm pitch angle directly. 
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by the author. 



